In thomasgaquiere/Maria: Simulate Different Forms of Logging on a Forest

INPUTS

inventory <- addtreedim(inventorycheckformat(Paracou6_2016), volumeparameters = ForestZoneVolumeParametersTable)

inventory <- treeselection(inventory, objective = 20, scenario ="manual", fuel = "2", diversification = TRUE, specieslax = FALSE,
                           objectivelax = FALSE, topography = DTMParacou, plotslope = PlotSlope, speciescriteria = SpeciesCriteria,
                           advancedloggingparameters = loggingparameters())$inventory


speciescriteria = SpeciesCriteria
scenario = "manual"
fuel = "2"
objective = 20
diversification = TRUE
directionalfelling = "2"
specieslax = FALSE
objectivelax = FALSE
topography = DTMParacou
advancedloggingparameters = loggingparameters()
# MainTrail

Succès-échec

# Selected <-  Accessible

successfail <- function(
  inventory,
  fuel = "2",
  directionalfelling = "2",
  advancedloggingparameters = loggingparameters()) {

  if (fuel == "0" && directionalfelling != "1"){
    inventory <- inventory %>% 
      mutate(TreeFellingOrientationSuccess = ifelse(Selected == "1", sample(c(1,0), size = 1, replace = F, prob = c(advancedloggingparameters$TreefallSuccessProportion, 1-advancedloggingparameters$TreefallSuccessProportion)), NA)) # Accessible = linked by 2ndtrails

    #Tree coordinates 
    if (any(inventory$TreeFellingOrientationSuccess == "1", na.rm = TRUE)) {

      TreefallSuccessCoord <- inventory %>% 
        filter(TreeFellingOrientationSuccess == "1") %>% 
        select(Xutm, Yutm)

      TreefallSuccess  <- st_multipoint(x = as.matrix(TreefallSuccessCoord))   # Create treefelling success
    }

    if (any(inventory$TreeFellingOrientationSuccess == "0", na.rm = TRUE)) {
      TreefallFailCoord <- inventory %>% 
        filter(TreeFellingOrientationSuccess == "0") %>% 
        select(Xutm, Yutm)

      TreefallFailure  <- st_multipoint(x = as.matrix(TreefallFailCoord))      # Create treefelling fails

    }

  }

  if (fuel =="1") {

    inventory <- inventory %>% 
      mutate(TreeFellingOrientationSuccess = ifelse(Selected == "1", sample(c(1,0), size = 1, replace = F, prob = c(advancedloggingparameters$TreefallSuccessProportion, 1-advancedloggingparameters$TreefallSuccessProportion)), NA)) # Selected = not yet linked by 2ndtrails, because 2ndtrails came after

    #Tree coordinates 
    if (any(inventory$TreeFellingOrientationSuccess == "1", na.rm = TRUE)) {

      TreefallSuccessCoord <- inventory %>% 
        filter(TreeFellingOrientationSuccess == "1") %>% 
        select(Xutm, Yutm)

      TreefallSuccess  <- st_multipoint(x = as.matrix(TreefallSuccessCoord))   # Create treefelling success
    }

    if (any(inventory$TreeFellingOrientationSuccess == "0", na.rm = TRUE)) {
      TreefallFailCoord <- inventory %>% 
        filter(TreeFellingOrientationSuccess == "0") %>% 
        select(Xutm, Yutm)

      TreefallFailure  <- st_multipoint(x = as.matrix(TreefallFailCoord))      # Create treefelling fails

    }
  }

  # Abattre les creux : 
  if (fuel =="2") {

    inventory <- inventory %>% 
      mutate(TreeFellingOrientationSuccess = ifelse(Selected == "1"| ProbedHollow == "1", sample(c(1,0), size = 1, replace = F, prob = c(advancedloggingparameters$TreefallSuccessProportion, 1-advancedloggingparameters$TreefallSuccessProportion)), NA)) # Selected = not yet linked by 2ndtrails, because 2ndtrails came after

    #Tree coordinates 
    if (any(inventory$TreeFellingOrientationSuccess == "1", na.rm = TRUE)) {

      TreefallSuccessCoord <- inventory %>% 
        filter(TreeFellingOrientationSuccess == "1") %>% 
        select(Xutm, Yutm)

      TreefallSuccess  <- st_multipoint(x = as.matrix(TreefallSuccessCoord))   # Create treefelling success
    }

    if (any(inventory$TreeFellingOrientationSuccess == "0", na.rm = TRUE)) {
      TreefallFailCoord <- inventory %>% 
        filter(TreeFellingOrientationSuccess == "0") %>% 
        select(Xutm, Yutm)

      TreefallFailure  <- st_multipoint(x = as.matrix(TreefallFailCoord))      # Create treefelling fails

    }
  }

  treefellingOutputs <- list(inventory = inventory,
                             TreefallSuccess = TreefallSuccess,
                             TreefallFailure = TreefallFailure)

  return(treefellingOutputs)

}

successfailOutputs <- successfail( #an inventory and 2 multipoints: 1 for success, 1 for fails
  inventory,
  fuel = "2",
  directionalfelling = "2",
  advancedloggingparameters = loggingparameters())

inventory <- successfailOutputs$inventory

felling <- function (inventory, 
                     SelectedTrees, 
                     fuel, 
                     directionalfelling, 
                     advancedloggingparameters = loggingparameters()
) {

  # Take selected trees
  # SelectedTrees <- inventory %>%
  #   filter(Selected == "1") %>%
  #   dplyr::select(idTree,DBH,TrunkHeight,TreeHeight,CrownHeight,CrownDiameter,Selected, Xutm, Yutm)

  # sp::coordinates(SelectedTrees) <- ~ Xutm + Yutm
  # 
  # sp::proj4string(SelectedTrees) <- raster::crs(topography)
  # 
  # SelectedTrees <- st_as_sf(as(SelectedTrees,"SpatialPoints"))

  # The Trunk
  pts = list(matrix(c(SelectedTrees$Xutm[1]-(SelectedTrees$DBH[1]/100)/2, SelectedTrees$Yutm[1], 
                      SelectedTrees$Xutm[1]-(SelectedTrees$DBH[1]/100)/2, SelectedTrees$Yutm[1] + SelectedTrees$TrunkHeight[1], 
                      SelectedTrees$Xutm[1]+(SelectedTrees$DBH[1]/100)/2, SelectedTrees$Yutm[1] + SelectedTrees$TrunkHeight[1], 
                      SelectedTrees$Xutm[1]+(SelectedTrees$DBH[1]/100)/2, SelectedTrees$Yutm[1],
                      SelectedTrees$Xutm[1]-(SelectedTrees$DBH[1]/100)/2, SelectedTrees$Yutm[1]) # the return
                    ,ncol=2, byrow=TRUE)) # DBH in cm to m

  (Trunk = st_polygon(pts))

  # The Crown
  dat <- data.frame(
    x = SelectedTrees$Xutm[1], #centroid location (x) same x than the tree
    y = SelectedTrees$Yutm[1] + SelectedTrees$TrunkHeight[1] + SelectedTrees$CrownHeight[1]/2, #centroid location (y) = trunk height + CrownHeight/2
    ex = SelectedTrees$CrownDiameter[1]/2, #Size along x-axis (CrownDiameter/2)
    ey = SelectedTrees$CrownHeight[1]/2, #Size along y-axis (CrownHeight/2)
    stringsAsFactors = FALSE
  )
  dat <- st_as_sf(dat, coords = c("x", "y")) # spacial dataframe

  Crown <- st_ellipse(pnt = dat, ex = dat$ex, ey = dat$ey) # Create the ellipse

  a <- st_difference(st_union(Crown, Trunk)) # -> multypolygon of which we keep only the points that do not overlap

  # Treefelling orientation

  # Create 2 polygons ScndTrail
  pol1 <- list(matrix(c(286503, 583134,
                        286503, 583240, 
                        286507, 583240, 
                        286507, 583134,
                        286503, 583134) # the return
                      ,ncol=2, byrow=TRUE)) 
  pol2 <- list(matrix(c(286650, 583134,
                        286650, 583240, 
                        286654, 583240, 
                        286654, 583134,
                        286650, 583134) # the return
                      ,ncol=2, byrow=TRUE))

  PolList = list(pol1,pol2) #list of lists of numeric matrices
  (ScndTrail = st_multipolygon(PolList))

  # Find the point (Trail) on the ScndTrail closest to the location of the tree (Foot)
  Foot <- st_point(c(SelectedTrees$Xutm[1],SelectedTrees$Yutm[1])) # tree foot point

  NearestPoints <- st_nearest_points(Foot, ScndTrail) # from the Foot of the tree to the ScndTrail (linestring)

  NearestPoint <- st_cast(NearestPoints, "POINT") # to have start (Foot) and end (Trail) points
  Trail <- NearestPoint[[2]] # the point (Trail) on the ScndTrail closest to the location of the tree (Foot)

  # Compute theta angle
  theta <- as.numeric(matlib::angle(c(Foot[1] - Foot[1], Foot[2] + SelectedTrees$TreeHeight[1]),
                                    c(Foot[1] - Trail[1], Foot[2] -Trail[2]), degree = TRUE))



  RotatePolygon <- function(p, angle, fixed) { # angle in the clockwise direction

    p.coords <- sf::st_coordinates(p)[,1:2] # Polygone coordinates extraction

    p.center <- suppressWarnings(sf::st_coordinates(fixed))

    rotate.coords <- function(xy, a, center) {

      co <- cos(-a * pi / 180)

      si <- sin(-a * pi / 180)

      adj <- matrix(rep(center, nrow(xy)), ncol=2, byrow=TRUE) # matrix with fixed point coordinates

      xy <- xy-adj

      cbind(co * xy[,1] - si * xy[,2],si * xy[,1] + co * xy[,2]) + adj
    }

    p.rotate <- rotate.coords(p.coords, a = angle, center = c(p.center[1], 
                                                              p.center[2]))

    Turned <- sf::st_sfc(sf::st_polygon(list(p.rotate))) # create the turned polygon

    return(Turned)

  }

  # For a random direction felling
  if(directionalfelling == "0" && (fuel !="1" |fuel !="2")){
    RandomAngle <- sample(c(0:180), size = 1)
    A <- st_difference(st_union(
      RotatePolygon(Trunk, angle = RandomAngle, fixed = Foot), # turned trunk
      RotatePolygon(Crown, angle = RandomAngle, fixed = Foot) # turned crown
    )) 
    B <- NULL
  }

  # To direct only to avoid damage to future and reserve trees
  # if (directionalfelling == "1"&& (fuel !="1" |fuel !="2")) {
  #   
  #   
  # }

  # To direct to avoid damage to future and reserve trees + trail orientation. Winching: Foot before.
  if(directionalfelling == "2"&& (fuel !="1" |fuel !="2")){
    A <- st_difference(st_union(
      RotatePolygon(Trunk, angle = 240 + theta, fixed = Foot), # turned trunk
      RotatePolygon(Crown, angle = 240 + theta, fixed = Foot) # turned crown
    ))
    B <- st_difference(st_union(
      RotatePolygon(Trunk, angle = 120 + theta, fixed = Foot), # turned trunk 
      RotatePolygon(Crown, angle = 120 + theta, fixed = Foot) # turned crown
    ))
  }

  # Fuel wood exploitation in the crowns
  if(fuel =="1" |fuel =="2"){

    TrailDist <- st_distance(Foot, Trail)

    # ADD SLOPE CRITERIA !!!

    if(TrailDist <= advancedloggingparameters$GrappleLength){ # <= 6m (= grapple length) -> winching by grapple
      A <- st_difference(st_union(
        RotatePolygon(Trunk, angle = theta + 60 , fixed = Foot), # turned trunk
        RotatePolygon(Crown, angle = theta + 60 , fixed = Foot) # turned crown
      ))
      B <- st_difference(st_union(
        RotatePolygon(Trunk, angle = 300 + theta, fixed = Foot), # turned trunk
        RotatePolygon(Crown, angle = 300 + theta, fixed = Foot) # turned crown
      ))
    } else { # > 6m -> winching by cable
      A <- st_difference(st_union(
        RotatePolygon(Trunk, angle = 240 + theta, fixed = Foot), # turned trunk
        RotatePolygon(Crown, angle = 240 + theta, fixed = Foot) # turned crown
      ))
      B <- st_difference(st_union(
        RotatePolygon(Trunk, angle = 120 + theta, fixed = Foot), # turned trunk 
        RotatePolygon(Crown, angle = 120 + theta, fixed = Foot) # turned crown
      ))
    }


  }

  FellingOuputs <- list(a = a,
                        Foot = Foot,
                        ScndTrail = ScndTrail,
                        NearestPoints = NearestPoints,
                        Trail = Trail,
                        A = A,
                        B = B)

  return(FellingOuputs)

}

Application

SelectedTrees <- inventory %>% 
  filter(Selected == "1") %>%
  dplyr::select(idTree,DBH,TrunkHeight,TreeHeight,CrownHeight,CrownDiameter,Selected, Xutm, Yutm)

EachTree <- list(SelectedTrees[1,],
                 SelectedTrees[2,],
                 SelectedTrees[3,])

AllTheFellingTrees <- lapply(EachTree, function(element)
  felling(inventory, element,  fuel = "2", directionalfelling = "2"))

a <- st_union(st_union(AllTheFellingTrees[[1]]$a, AllTheFellingTrees[[2]]$a), AllTheFellingTrees[[3]]$a)

Foot <- st_union(st_union(AllTheFellingTrees[[1]]$Foot, AllTheFellingTrees[[2]]$Foot), AllTheFellingTrees[[3]]$Foot)

ScndTrail <- st_union(st_union(AllTheFellingTrees[[1]]$ScndTrail, AllTheFellingTrees[[2]]$ScndTrail), AllTheFellingTrees[[3]]$ScndTrail)

NearestPoints <-  st_union(st_union(AllTheFellingTrees[[1]]$NearestPoints, AllTheFellingTrees[[2]]$NearestPoints), AllTheFellingTrees[[3]]$NearestPoints)
# Trail <- Trail
A <- st_union(st_union(AllTheFellingTrees[[1]]$A, AllTheFellingTrees[[2]]$A), AllTheFellingTrees[[3]]$A)

B <- st_union(st_union(AllTheFellingTrees[[1]]$B, AllTheFellingTrees[[2]]$B), AllTheFellingTrees[[3]]$B)


ggplot() + # plot 2 polygones
  geom_sf(data = a, colour = "red") + # default tree position (polygone)
  geom_sf(data = Foot) + # the foot of the tree (point)
  geom_sf(data = ScndTrail) + # the 2nd trail (polygone)
  geom_sf(data = NearestPoints) + # the shortest way from the foot to the trail
  # geom_sf(data = Trail) + # the point on the ScndTrail closest to the location of the tree (Foot)
  geom_sf(data = A, colour = "green") + # Foot before
  geom_sf(data = B, colour = "green")